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Abstract 

We provide a mathematical analysis and a numerical framework for mag¬ 
netoacoustic tomography with magnetic induction. The imaging problem is to 
reconstruct the conductivity distribution of biological tissue from measurements 
of the Lorentz force induced tissue vibration. We begin with reconstructing from 
the acoustic measurements the divergence of the Lorentz force, which is acting 
as the source term in the acoustic wave equation. Then we recover the electric 
current density from the divergence of the Lorentz force. To solve the nonlin¬ 
ear inverse conductivity problem, we introduce an optimal control method for 
reconstructing the conductivity from the electric current density. We prove its 
convergence and stability. We also present a point fixed approach and prove 
its convergence to the true solution. A new direct reconstruction scheme in¬ 
volving a partial differential equation is then proposed based on viscosity-type 
regularization to a transport equation satisfied by the electric current density 
field. We prove that solving such an equation yields the true conductivity dis¬ 
tribution as the regularization parameter approaches zero. Finally, we test the 
three schemes numerically in the presence of measurement noise, quantify their 
stability and resolution, and compare their performance. 


Mathematics Subject Classification (MSC2000): 35R30, 35B30. 

Keywords: electrical conductivity imaging, hybrid imaging, Lorentz force induced tissue vibration, 
optimal control, convergence, point fixed algorithm, orthogonal field method, viscosity-type regu¬ 
larization. 


1 Introduction 

The Lorentz force plays a key role in magneto-acoustic tomographic techniques 
[36]. Several approaches have been developed with the aim of providing electri- 
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cal impedance information at a spatial resolution on the scale of ultrasound wave¬ 
lengths. These include ultrasonically-induced Lorentz force imaging mm and 
magneto-acoustic tomography with magnetic induction [laEZ]. 

Electrical conductivity varies widely among soft tissue types and pathological 
states [231 [35] and its measurement can provide information about the physiological 
and pathological conditions of tissue [HI . Acousto-magnetic tomographic techniques 
have the potential to detect small conductivity inhomogeneities, enabling them to 
diagnose pathologies such as cancer by detecting tumorous tissues when other con¬ 
ductivity imaging techniques fail to do so. 

In magnetoacoustic imaging with magnetic induction, magnetic fields are used 
to induce currents in the tissue. Ultrasound is generated by placing the tissue in 
a dynamic and static magnetic field. The dynamic field induces eddy currents and 
the static field leads to generation of acoustic vibration from Lorentz force on the 
induced currents. The divergence of the Lorentz force acts as acoustic source of 
propagating ultrasound waves that can be sensed by ultrasonic transducers placed 
around the tissue. The imaging problem is to obtain the conductivity distribution 
of the tissue from the acoustic source map; see [SllEaiMlISlIll- 

This paper provides a mathematical and numerical framework for magnetoacous¬ 
tic imaging with magnetic induction. We develop efficient methods for reconstruct¬ 
ing the conductivity in the medium from the Lorentz force induced vibration. For 
doing so, we first estimate the electric current density in the tissue. Then we design 
efficient algorithms for reconstructing the heterogeneous conductivity map from the 
electric current density with the ultrasonic resolution. 

The paper is organized as follows. We start by describing the forward problem. 
Then we reconstruct from the acoustic measurements the divergence of the Lorentz 
force, which is acting as the source term in the acoustic wave equation. We recover 
the electric current density from the divergence of the Lorentz force, which reduces 
the problem to imaging the conductivity from the internal electric current density. 
We introduce three reconstruction schemes for solving the conductivity imaging 
problem from the internal electric current density. The first is an optimal control 
method. One of the contributions of this paper is the proof of convergence and 
stability of the optimal control approach provided that two magnetic excitations 
leading to nonparallel current densities are employed. Then we present a point fixed 
approach and prove that it converges to the true conductivity image. Finally, we 
propose an alternative to these iterative schemes via the use of a transport equation 
satisfied by the internal electric current density. Our third algorithm is direct and 
can be viewed as a PDF-based reconstruction scheme. We test numerically the 
three proposed schemes in the presence of measurement noise, and also quantify 
their stability and resolution. 

The feasibility of imaging of Lorentz-force-induced motion in conductive sam¬ 
ples was shown in [21j . The magnetoacoustic tomography with magnetic induction 
investigated here was experimentally tested in [33l IMj , and was reported to produce 
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conductivity images of quality comparable to that of ultrasound images taken under 
similar conditions. Other emerging hybrid techniques for conductivity imaging have 
also been reported in [DElEllillllliaiMlIMlIlS]. 

2 Forward problem description 

2.1 Time scales involved 

The forward problem in magnetoacoustic tomography with magnetic induction (MAT- 
MI) is multiscale in nature. The different phenomena involved in the experiment 
evolve on very different time scales. Precisely, there are three typical times that 
appear in the mathematical model for MAT-MI. 

• The first one is the time needed for an electromagnetic wave to propagate in 
the medium and is denoted by Tem- Typically, if the medium has a diameter 
of 1cm, we have Tem ~ 

• The second characteristic time length, denoted by Tpuise of the experiment is 

the time width of the magnetic pulse sent into the medium. Since, the time- 
varying magnetic field is generated by discharging a capacitor, Tpuise is in fact 
the time needed to discharge the capacitor such that Tpuise ~ [33]. 

• The third characteristic time, Tsound, is the time consumed by the acoustic wave 
to propagate through the medium. The speed of sound is about 1.5 • 10^m.s“^ 
so Tsound ~ 6//S for a medium of 1cm diameter. 

2.2 Electromagnetic model 

Let (ej)i=i^ 2,3 be an orthonormal basis of M^. Let be a three-dimensional bounded 
convex domain. The medium is assumed to be non magnetic, and its conductivity 
is given by a (the question of the regularity of a will arise later). Assume that the 
medium Q is placed in a uniform, static magnetic field in the transverse direction 
Bq = ^063. 


2.2.1 The magnetoquasistatic regime 

At time t = 0 a second time varying magnetic field is applied in the medium. The 
time varying magnetic field has the form Bi(x,t) = Bi{x)u{t)e 3 . Bi is assumed 
to be a known smooth function and u is the shape of the stimulating pulse. The 
typical width of the pulse is about l//s so we are in presence of a slowly varying 
magnetic-field. This regime can be described by the magnetoquasistatic equations 
|28| . where the propagation of the electrical currents is considered as instantaneous, 
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but, the induction effects are not neglected. These governing equations in x 
are 

V-B = 0, (2.1) 

9B 

VxE = - —, (2.2) 

and 

V-J = 0, (2.3) 

where eo is the permittivity of free space, B is the total magnetic field in the medium 
and E is the total electric field in the medium. Ohm’s law is valid and is expressed 
as 

J = crE in n X M_|_, (2.4) 

where cr is the electrical conductivity of the medium. For now on, we assume that 
a G L“^(n), where 

L^hi^) := {/ G : a < / < 6 in O', / = uq in 0 \ IT} 

with (To, a, and b being three given positive constants, 0 < a < b, and 12'd 12. 

As in [28], we use the Coulomb gauge (V • A = 0) to express the potential 
representation of the fields B and E. The magnetic field B is written as 


B = V X A, 


(2.5) 


and the electric field E is then of the form 

dA 


E = -VE- 


dt 


in 12 X 




( 2 . 6 ) 


where V is the electric potential. Writing A as follows: 

A{x,t) = Ao(x) + Ai{x)u{t), 


where Aq and Ai are assumed to be smooth. In view of (2.3) and (2.6), we look for 
V{x,t) of the form V{x,t) = V{x)u'{t) with V satisfying 


V • aW = —V • (tAi in 12 X M. 


The boundary condition on V can be set as a Neumann boundary condition. Since 
the medium 12 is usually embedded in a non-conductive medium (air), no currents 
leave the medium, i.e., J • = 0 on 512, where i/ is the outward normal at 512. To 

make sure that the boundary-value problem satisfied by V is well posed, we add the 
condition /qV = 0. We have the following boundary value problem for V: 


V ■ aVV = - V • o-Ai in 12, 
dV 

a-;— = — (tAi • IV on 512, 
ou 


V =0. 


/n 


(2.7) 
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2.3 The acoustic problem 
2.3.1 Elasticity formulation 

The eddy currents induced in the medium, combined with the magnetic field, create 
a Lorentz force based stress in the medium. The Lorentz force f is determined as 


f = JxB infix M+. 


( 2 . 8 ) 


Since the duration and the amplitude of the stimulation are both small, we 
assume that we can use the linear elasticity model. The displacements inside the 
medium can be described by the initial boundary-value problem for the Lame system 
of equations 


pd‘^vL 




VAV • u — V • /rV^u = J X B 


u(x, 0) 


(9u 


d\i 

dn 

(x,0) = 0 


in fl X M_|_, 
on dQ X M+, 

in fl. 


(2.9) 


where A and p are the Lame coefficients, p is the density of the medium at rest, 
and V*u = (Vu + V^u)/2 with the superscript T being the transpose. Here, djdn 
denotes the co-normal derivative defined by 

dw. 

— = A(V • u)^/ -|- 2uV^ui>' on 5H, 
on 

where v is the outward normal at dVt. The functions A, p, and p are assumed to be 
positive, smooth functions on H. 

The Neumann boundary condition, du/dn = 0 on Sfl, comes from the fact that 
the sample is embedded in air and can move freely at the boundary. 


2.3.2 The acoustic wave 


Under some physical assumptions, the Lame system of equations (2.9) can be re¬ 
duced to an acoustic wave equation. For doing so, we neglect the shear effects in 
the medium by taking p = 0. The acoustic approximation says that the dominant 
wave type is a compressional wave. Equation (2.9) becomes 

pd^u — VAV •u = JxB infix IR+, 


clu 

— = 0 on 5fl X 
on 




( 2 . 10 ) 


5u 


Introduce the pressure 


u(x,0) = —(x,0) = 0 in fl. 


p = AV • u in fl X 
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Taking the divergence of (2.10) yields the acoustic wave equation 


1 1 1 

— —w — V • -Vj? = V • - (J X B) in n X M+, 

A p p 

p = 0 on dVL X IR+, 
dv 

p(x,0) =-^(x,0) = 0 in 


( 2 . 11 ) 


We assume that the duration Tpuigg of the electrical pulse sent into the medium 
is short enough so that p is the solution to 


1 d'^v 1 

= f{x)6t=o in n X M+, 

A ot^ p 

p = 0 on do, X M+, 
p(a:,0) = ^(x,0) = 0 in 


( 2 . 12 ) 


where 


/'-* pulse 

f(x)= / V ■ (-J(x,t) X'B(x,t))dt. 
Jo P 


(2.13) 


Note that acoustic wave reflection in soft tissue by an interface with air can be 
modeled well by a homogeneous Dirichlet boundary condition; see, for instance, (4lj . 
Let 

dv 

g(x,t) = —{x,t), V(x,t) G cin X M+. 
dv 

In the next section, we aim at reconstructing the source term / from the data g. 

3 Reconstruction of the acoustic source 

In this subsection, we assume that A = Aq + 5A and p = po + 5p, where the functions 
6X and 6p are such that ||4A|<C Aq and ||4p||Loo(r2) po- We assume that 

— 


A, Ao, p, and po ai'e known and denote by cq = y ^ the background acoustic speed. 

Based on the Born approximation, we image the source term /. To do so, we first 

consider the time-harmonic regime and define r‘^ to be the outgoing fundamental 
2 

solution to A -I- 

Cfl 


OJ 


Ax + ^ r^(a;,2/) = 6y{x), 


(3.1) 


subject to the Sommerfeld radiation condition: 


1 

\x \2 




0 , 


oo. 
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We need the following integral operator : L?‘{d^) —?■ L?‘{d^) given by 

w^(a:,2/)0(y)ds(2/). 

Jan di'ix) 

2 

Let Gq be the Dirichlet Green function for A + ^ in fl, i.e., for each y £ Q, Gq{x, y) 
is the solution to 


(^^x +Gf^{x,y) = 5y{x), x£n, 
G'^{x,y) = 0, x £ do.. 


Let p denote the Fourier transform of the pressure p and g the Fourier transform 
of g. The function p is the solution to the Helmholtz equation: 

^2 

p(x,a;) + V • —-Vp(x,a;) =/(x), x G H, 

A(x) p{x) 

p{x, uj) = 0, X £ dQ. 


Note that / is a real-valued function. 

The Lippmann-Schwinger representation formula shows that 

P{x,uj) = Jj,^^-l)Vp{y,u})-VGfi{x,y)dy-uj^ ^)p{y,uj)Gfi{x,y)dy 

+Po [ f{y)Gn{x,y)dy. 

Jn 

Using the Born approximation, we obtain 

p{x,uj) PS —- f 6p{y)Vpo{y,uj) ■VGfi{x,y)dy+ f ^^^po{y,u;)G‘i{{x,y) dy 

Po Jn Cq Jq Ao 

+Po f f{y)Gn{x,y)dy 
Jn 

for X G H, where 


Poix,uj):=po fiy)Gfi{x,y)dy, x G H. 
Jn 


Therefore, from the identity [HI Eq. (11.20)] 


( 2 ^ + ra^ 

it follows that 


801 

dv. 


{■,y) 


(x) = 


9F^ 

dv{x) 


(x, y), X G 9H, y G H, 


1 

po Jn 


8p{y)VpQ{y,u) ■ V 


9F‘^(x,y) f 6X{y) 


dv{x) 


dy + 


Cq Jn ^0 


po(y,w) 


5F‘^(x, y) 
dv{x) 


dy 


, f ^dT-{x,y) ^ 
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for X G do.. 
Introduce 


I{z,uj) : = 


Ion . 


r"(i,2)(J/+(CS)*)|9](i,w)-r“(i.3)(l; + (CS)*)|9l(i^,^) 


ds{x) 


for z € Q. 

We recall the Helmholtz-Kirchhoff identity m Lemma 2.32] 


Ian 


—, .dT‘^(x,y) .dT‘^(x,y) 

r^(x,z) ^ -r‘^(x,z) ^ 


dv{x) 


dv{x) 


ds{x) = 2i3‘mr‘^(2;, y). 


We also recall that / is real-valued and write / -I- 5f. Given I{z,uj) we solve 

the deconvolution problem 


2ipQ QmT‘^{z,y)f'^'^\y)dy = I{z,uj), z € n, 


(3.2) 


in order to reconstruct with a resolution limit determined by the Ra yleig h 
criteria. Once is determined, we solve the second deconvolution problem (3.3) 


2ipo QmT‘^{z,y)6f{y)dy = 6I{z,io), z € n, 
Jn 

to find the correction 6f. Here, 


(3.3) 


61{z,uj) := 


Ian . 


r‘^{x, z)6g{x,uj) —r‘^{x,z)6g{x,uj) 


ds{x) 


with 


1 /■ f dX{y) ^(n\ dr‘^(x,y) 

5g{x,u}) = — / 6p{y)Vp^°\y,uj)-V——^dy+^ / -—^p^°\y,u}) . dy, 

Po Jn dv{x) Vq Ao dv{x) 


and 


p^^\x,uj) := po f^^\y)GQ{x,y)dy, x G H. 
Jn 


Since by Fourier transform, g is known for all oj G M+, I{z,uj) can be computed for 
all uj G M+. Then from the identity |TT1 Eq. (1.35)] 


2 

TT 


/ a;3‘mr^(x, z) doj = —5z{x), 

J M-|- 

where 5z is the Dirac mass at z, it follows that 

f^^\z) =-. - [ u!l{z,oj)duj and 6 f{z) = - - [ uj5I{z,uj) duj. 

Jm+ i-t^Po Jr+ 


We refer to |19l [20] and the references therein for source reconstruction approaches 
with finite set of frequencies. 


















4 Reconstruction of the conductivity 


We assume that we have reconstructed the pressure source / given by (2.13). We also 
assnme that the sample O is thin and hence can be assimilated to a two dimensional 
domain. Further, we suppose that 0 C vect ( 61 , 02 ). Here, vect ( 61 , 62 ) denotes the 
vector space spanned by 61 and 62 . Recall that the magnetic fields Bq and Bi are 
parallel to 63 . We write J{x,t) = J(x)tt'(t). In order to recover the conductivity 
distribution, we start by reconstructing the vector field J(x) in 0. 


4.1 Reconstruction of the electric current density 

4.1.1 Helmholtz decomposition 

Let H^{Q) := {n G L^(n) : Vn G L^(fl)}. Let llg(n) be the set of functions in 
with trace zero on dQ and let be the dual of Hq{Q,). 

We need the following two classical results. 

Lemma 4.1. If a G L“(,(fl) then the solution V of 
hence, the electric current density J belongs to if {11). 

The following Helmholtz decomposition in two dimensions holds [lOj . 

Lemma 4.2. //f is a vector field in L‘^{Q), then there exist two functions v G H^{Q) 
and w G H^{Il) such that 

f = Vn + curl w. (4-1) 

The differential operator curl is defined by curl w = {—d 2 W,diw). Furthermore, if 
V • f G if {Ft), then the potential v is a solution to 


'2.1) belongs to H^{Fl) and 


—Av = V • f in Ft, 
dv 

— = 1-1/ on oil, 
ov 


(4.2) 


and w is the unique solution of 


curl w ■ curl = / (f — Vn) • curl Mcf G H{FL), 

Jo, 


(4.3) 


where H{Fl) = {<)) G L^(fl), V x ^ L^(fl), V ■ = 0}. The problem can be written 

in strong form in H~^{Fl): 


—Aw = curl f 
tc = 0 


in n, 
on dFl, 


where the operator curl is defined on vector fields by curl f = —d 2 fi + dif 2 . 
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We apply the Helmholtz decomposition (4.1) to the vector field J G and 

get the following proposition. 


Proposition 4.3. There exists a function w € H such that 

J = curl w, 


(4.4) 


and w is the unique solution of 


Recall (2.3): 



curl J 
w = 0 


in 

on dQ. 


V-J = 0, 


together with the fact that no current leaves the medium 


J • iz = 0 on 0H. 


(4.5) 


Since u is a solution to (2.7), v has to be constant. So, in order to reconstruct J one 
just needs to reconstruct w. 


4.1.2 Recovery of J 

Under the assumption |Bi| <C |Bo| in x M+ and \5p\ <C po in H, the pressure 
source term / defined by (2.13) can be approximated as follows: 

f(x) ^ Iv • (J(x) X Bo)(u(rp,i,e) - u(0)), 

po 

where we have used that J(x,t) = i{x)u'{t). 

Since Bq is constant we get 


V • (J(x) X Bq) = (V X j) • Bq = |Bo|curl J. 


Now, since Bq is known, we can compute w as the unique solution of 

_Po/_ ^ 

Bo|(^^(rpuise) - u(0)) ™ ’ (4.6) 

rc = 0 on dTl, 

and then, by Proposition |4.3[ compute J by J = curl w. 

Note that since the problem is reduced to the two dimensional case, J is then 
contained in the plane Bq with T denoting the orthogonal. 
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4.2 Recovery of the conductivity from internal electric current den¬ 
sity 

In this subsection we denote by ct* the true conductivity of the medium, and we 
assume that <t* G with 0 < a < b, i.e., it is bounded from below and above 

by positive known constants and is equal to some given positive constant do in a 
neighborhood of dQ. 


4.2.1 Optimal control method 


Recall that Ai is defined by V • Ai = 0, Bi{x)e^ = V x Ai(x). Define the following 
operator T: 

a I—;■ 

with 

( V • aVU = — V • (tAi in D, 


F[a] := U { 


The following lemma holds. 


A an 

= — dAi • u on oil, 

ov 


(4.7) 


t/=0. 




Lemma 4.4. The operator T is Frechet differentiable. For any a G and h 

such that a + h G we have 


dT[a]{h) :=q{ 


rv • dVg = — V • /lAi — V • hVF[a] in D, 

=0 on dD, 

on 

[ q=0. 


(4.8) 


Proof. Denote by r the function F\a + /i] — F[(t] — q. The function r belongs to 
and satisfies the following equation in D: 

V • dVr = V • /iV {F[a] - F[a + h ]), 

together with the boundary condition 

dr 


du 


= 0 on dD, 


and the zero mean condition = 0. We have the following estimate: 


< -||/r|Uoo(n)||V {F[a] - F[a + h]) 
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Since F[a] — F\a + h\ satisfies 


V • (cjV {F[a] - T\a + h])) = -V • {hVF[a + h])+V- (hAi) 
with the boundary condition 

^{F[a + h]-F[a]) = 0, 

and the zero mean condition {F[a + h]— J^[cj]) = 0. We can also estimate the 
L^-norm of V {F[a + h]— F[a]) as follows; 

||V {F[a + h]— F[a]) \\L'2{n) < “ll^llL“(n) + ^]llL2(f7) + ||Ai||2,2(q)) . 

Therefore, we can bound the ff^-norm of F{a + h\ independently of a and h for 
11/i| I LOO small enough. There exists a constant C, depending only on fl, a, b, and 
Ai, such that 

||V.T[c7 + /i]||l2(^)<C. 

Hence, we get 

\\ViF[a + h]-F[a]) ||l2(o) < Iwhh^in) (C + ||Ai||L2(n)) , 
and therefore, 

l|Vr||L2(j^) < C||/i|||oo(f^), 

which shows the Frechet differentiability of F. □ 

Now, we introduce the misfit functional: 

a ^ J[a] = l [ \a (V.TH + Ai) - Jp, 

Jn 

Lemma 4.5. The misfit functional J is Frechet-differentiable. For any a £ 
we have 


dj[a] = (cjV.T[cj] + uAi - J) • {VF[a] + Ai) + Vs • (Ai + V.T[cj]), 

where s is defined as the solution to the adjoint problem: 

V • aVs =V • (cj^VT'[cj] + ci^Ai — cjj) in H, 
ds 

a— =0 on dil, 

ou 

f s=0. 


(4.10) 
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Proof. Since J- is Frechet-difFerentiable, so is ff. For any a G and h such 

that a + h G we have 


dj[a]{h)= / {aVT[a]+aAi-3)-{hV{J^[a] + Ai) + aV{dT[a]{h))) 
Jn 


Multiplying (4.10) by dP'[a]{h) we get 


/ (t(cjV-F[u]+ (tAi - J) • VdJ^[fj](/i) = / aVs-VdT[a]{h). 
In Jn 


On the other hand, multiplying (4.8) by s we obtain 


[ aVs-VdJ^[a]ih) = [ Ws • (Ai + . 

Jn Jn 


So we have 


dj[a]{h)= / h (cjV-F[fT] + (TAi-J)-(VJ^[fj]+Ai)+Vs-(Ai + VJ'[(7]) 

Jn i 

and the proof is complete. 


□ 


Lemma 4.5 allows us to apply the gradient descent method in order to minimize 
the discrepancy functional J'. Let (T(o) be an initial guess. We compute the iterates 

cT(n+i) = r[cr(„)] - fidJ[T[a^n)]], Vn G N, (4.11) 

where fi > 0 is the step size and T[f] = min{max{/, o}, b}. 

In the sequel, we prove the convergence of ( 4.11| ) with two excitations. Let 
and correspond to two different excitations a[^^ and A^^\ Assume that 
j{i) X 7 ^ 0 in 0. Let : a i—)• aV -|- A^*^^ — Jj, where is dehned 

I I J I 

by (4.8) with Ai = A) ' for i = 1,2. The optimal control algorithm (4.11) with two 
excitations is equivalent to the following Landweber scheme given by 


= r[cj(„)] - /rda*[^[T[o-(„)]]], Vn G N, (4.12) 


where G[a] = (^(^)[cr], ^(^^[cj])'^. 

Following |15j . we prove the convergence and stability of (4.12) provided that 


two magnetic excitations leading to nonparallel current densities are employed. 


Proposition 4.6. Let and correspond to two different excitations. Assume 
that X 7^ 0 in n. Then there exists rj > 0 such that if ||cr(o) — < V; 

then I|cr(„) — <7*|—)■ 0 as n —)• +oo. 
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Proof. According to m, it suffices to prove that there exists a positive constant C 
such that 

\\dG[a]{h)\\Hii^)>C\\h\\j,r^^) (4.13) 

for all h G such that u + /i G L“^(Q). We have 

dg^^[a]{h) = /ijW + aVdF^^[a]{h). 


Therefore, 

V • dg«[cj](h) = 0, dg^^[a]{h) • 1 / = 0, 

and 

V X (-dg«[fT](/l)) = W X (-J«) +f7V/l X J«. 
G a 

Since V x x 63 = 0 and x / 0, it follows that 

2 

i=l 


which completes the proof. 


□ 


Let F[a\ = 

positive constant C such that 


Note that analogously to (4.13) there exists a 


\\dF[a\{h)\\H^^^>C\\h\\Hi(n) 

for all h G Hq{Q) such that a + h € L“^(n), provided that x ^ 0 in 11. 


4.2.2 Fixed point method 

In this subsection, we denote by ct* the true conductivity inside the domain H. We 
also make the following assumptions: 

• 3c > 0, such that |Bi| > c in 11; 

• cr G C°’“(n), a g]0, 1[; 

• cj* = uq in an open neighborhood of dll. 

From the unique continuation principle, the following lemma holds. 

Lemma 4.7. The set {x G H, J(x) = 0} is nowhere dense. 
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The interior data is J = a* [VT'[<t*] + Ai], One can only hope to recover cr* at 
the points where J 7 ^ 0. Even then, we can expect any type of reconstruction to 
be numerically unstable in sets where J is very small. Assume that J is continuous 
and let e > 0 and xq be such that |J(xo)| > 2 e. We define Og to be a neighborhood 
of xo such that for any x G Og, |J(x)| > e. One can assume that Og is a domain 
without loosing generality. Now, introduce the operator Tg as follows: 

a I— >Fe[(T] := Eg, 

where I 4 satisfies the following equation: 


V • crVEg = — V • ((tAi) in Og 

(T —— = — crAi ■ u + J ■ u on asZg 
ov 


(4.14) 


K=0, 


/He 


where v denotes the outward normal to clOg. Note that J • = 0 since V • J = 0 

in Og. 

We also define the nonlinear operator by 




cr I— >QeW] '■= cr-^ 


(cjVEg[cj] +uAi) • J 


(4.15) 


Lemma 4.8. The restriction of cr* on fig is a fixed point for the operator Q^. 

Proof. For the existence it suffices to prove that T'g [cT*|nE] = • Denote by 

14 = T'[cr*]. We can see that E* satisfies 


V • c 7*VE* = —V • (uAi) in Dg. 
Taking the normal derivative along the boundary of Dg, we get 

a —— = —(tAi ■ v' + j ■ 1 / on aSZg. 
ou 


From the well posedness of (4.14), it follows that 


So, we arrive at 


E*|^^ = J4k*|^J+c, cGM. 

Ge ~ *^*10 ■ 


□ 
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We need the following lemma. We refer to m for its proof. 

Lemma 4.9. Let 0 C be a bounded domain with Lipschitz boundary. For each 
g G there exists at least one v G with = g in the sense of the 

distributions and 

l|v||L2(f^) < C\\g\\H-^Q) 
with the constant C depending only on fl. 

The following result holds. 

Lemma 4.10. If \\Ai\\ 1 ^ 2 is small enough, then the operator is a contraction. 
Proof. Take cJi and a 2 in We have 

\ge[ai]{x) - a£[o- 2 ](x)| = 

X \{al{x)VVe[ai]{x) - a 2 {x)VVs[a 2 ]{x) + (o-f(x) - cr^ix)) Ai(x)) • J(x)| , 
which gives, using the Cauchy-Schwartz inequality: 

\SeWl]ix) - ge[o'2]{x)\ < ^ 

X I (al{x)VVe[ai]{x) - a 2 {x)We[a 2 ]{x) + {al{x) - al{x)) Ai(x)) | . 

The right-hand side can be rewritten using the fact that |o'j(x)| < b for i = 1,2, and 
hence, 


\Qe[(Tl]{x) - ge[o-2](x)\ < ^ 

X [|cTi(a:)V14[(Ti](x) - <T2{x)We[a2]{x)\ + |((Ti(x) - CT2 (x)) Ai(x)|] . (4.16) 
Now, consider the function v = (TiV14[cji] — fT 2 V 14 [(T 2 ]. We get 


V • V = —V • [(fJi — CJ 2 ) Ai] in dLl^, 


along with the boundary condition v • iz = 0 on 
a constant C depending only on fig such that 


Using Lemma 4.9 there exists 


I|v||l2(q^) < CIIV • [((Ti - CJ 2 ) Ai] ||_H--i(ng), 

which shows that 

l|v||L2(t2^) < C\\ {ai - 0 - 2 ) Ai ||L2(f2^). 
Using Cauchy-Schwartz inequality: 


v|| 2 , 2 (q^) < C||cri - 0 - 2 ||i 2 (f 2 ^)||Ai ||L 2 (f 2 ^). (4.17) 
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Putting together (4.16) with (4.17), we arrive at 


||^£[cri] - GeW2]\\L2(n^) < {C + 1) ^ || Ai ||i2(Q^) ||(Ti - Cr2\\L2{n^)- 
The proof is then complete. 

□ 


The following proposition shows the convergence of the hxed point reconstruction 
algorithm. 

Proposition 4.11. Let cr(^) G (L^(n£))^ be the sequence defined by 


<^(o) — 1 ) 

cr(„+i) = max (min {GeW{n)],b) , a) , Vn G N. 


(4.18) 


If II Ai is small enough, then the sequence is well defined and ai^^) converges 
to cr*|^^ in L‘^{Lls)- 

Proof. Let {X,d) = ^L“^(n£), || • ||L 2 (f^^)^. Then, {X,d) is a complete, non empty 
metric space. Let Tg be the map defined by 


[<t] := max (min [u] ,b), a) 


Using Lemma 
small enough 


4.10 


we get that 


Tg is a contraction. 


provided that ||Ai|| 2 , 2 (q^) is 


We already have the existence of a hxed point given by Lemma |4.8t 
and therefore, Banach’s hxed point theorem gives the convergence of the sequence 
for the norm over and the uniqueness of the hxed point. □ 


4.2.3 Orthogonal field method 


In this section we present a non-iterative method to reconstruct the electrical con¬ 
ductivity from the electric current density. This direct method was hrst introduced 
in m and works with piecewise regularity for the true conductivity a* in the case 
of a Lorentz force electrical impedance tomography experiment. However, the prac¬ 
tical conditions are a bit different here and we have to modify the method to make 
it work in the present case. 

We assume in this section that cr* G C'^’“(n), a g]0, 1]. The helds J = (Ji, J 2 ) 
an d Ai are assumed to be known in U. Our goal is to reconstruct 14 the solution 
of (2.7) in Then, computing for |J| nonzero will give us Recall 


that J = curl w where w is dehned by equation (4.6) 


Definition 4.1. We say that the data f on the right hand side of {4-6) is admissible 
if f > 0 or f < 0 in Q and if the critical points of w are isolated. 
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Introduce F = (—J 2 , Ji)^ the rotation of J by It is worth noticing that the 
true electrical potential I* is a solution of 


r F • VK = -F • Ai 
1^=0 
[ 14 = 0. 

Jn 


in fl, 
on dQ, 


(4.19) 


Equation (4.19) has a unique solution in and this solution is the true potential 

F*. 

The following uniqueness result holds. 

Proposition 4.12. If U £ H^{Q) is a solution of 

in Q, 


(F-VU = 0 

on 

[ U = 0, 


on dQ, 


(4.20) 


then U = D in Ft. 

Proof. We use the characteristic method (see, for instance, 
For any xq G Ft, consider the Cauchy problem: 


dX 


= F{X{t)), t£ 


) for solving (4.20). 


(4.21) 


dt 

A(0) = xo G n. 

We call the set {x{t),t G M} the integral curve at xq. Since a G C‘’’"(n), F G C^’"(n). 
Then, we can apply the Cauchy-Lipschitz theorem and get global existence and 
uniqueness of a solution to (4.21). Denote by T the upper bound on the domain 
size of integral curves. Now, assume that U G H^{Fl) is a solution of (4.20). Since 
J = curl w, F can be written as 

F = —Vw in Ft. 


Equation (4.21) reduces to the following gradient flow problem: 

^ = -Vw{Xit)), t£R, 

X(0) = xo G D. 


(4.22) 


Using [35], we know that there are hnitely many isolated critical points pi,... ,pn, 
for w on Q. It is also known (see [Ml p. 204]) that since the sets w~^ (] — 00 , cj) are 
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compact for every c G M, limt_^oo exists and is equal to one of the equilibrium 
points pi,... ,pn- Now, for every i, we define Qi the set of points xq G such that 
the solution of (4.22) converges to pi. Therefore, we have 12 = 


Now, for any i consider xq G 12*, and X G ([0,r[, 12) the solution of (4.22). 
We define f G (]R''',M) by f{t) = U{X{t)). The function / is differentiable on M+ 
and f'{t) ='VU{X{t)) ■ F{X{t)) = 0. Hence, / is constant. We have 

U{xo) = /(O) = lim f{t) = U{pi) = Ci G M. 

t^OO 

So, U is constant equal to c* in 12j. The regularity of U implies that Vz,j G [[l,u]], 
Cj = Cj. Therefore U is constant on 12 and the zero integral condition yields 

U = 0 in 12. 


This shows the uniqueness of a solution to (4.20) and thus, concludes the proof. □ 


In order to solve numerically (4.19), we use a method of vanishing viscosity m- 
The field Ai is known and we can solve uniquely the following problem: 


V 


{rjl + FF^) 

dU^v) 


dv 

= 0 , 


=-V•FF^Ai 


= — Ai • u 


in 12, 
on 912, 


(4.23) 


for some small rj > 0. Here, I denotes the 2x2 identity matrix. 


Proposition 4.13. Let cj* be the true c onductivity. Let 14 be the solution to (^-D 
with a = cr*. The solution Lf^'^'i of (4-23} converges strongly to V* in H^{Ll) when rj 
goes to zero. 

Proof. We can easily see that — 14 is the solution to 

V ■ ,, 

artlv) 

(4.24) 


= -r?AI4 

in 12, 

o 

II 

CO 

on 912, 

/ = 0. 

Jn 



Multiplying (4.24) by and integrating by parts over 12, we find that 

p [ [ |F-=r/ [ + [ U^^^Ai-u, (4.25) 

Ju Jci J Cl JdCl 
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since —-— = 0 and —— = —Ai • i/. Therefore, we have 
ou ou 

l|i/‘''>ll?,.(n) < l|£''’”llH.(nllKllH.(n) +C||cf’'ll„.,n), 

where C depends only on and Ai. This shows that the sequence {U^'^^)n>o is 
bounded in Using Banach-Alaoglu’s theorem we can extract a subsequence 

which converges weakly to some u* in We multiply (4.24) by u* and integrate 

by parts over Q, to obtain 


in 


[ (F • (F • Vu*) =r] 

Jn 

Taking the limit when r] goes to zero yields 

I|F-Vu*||i2(^)=0 


VU* -Vu* - / • Vm* + / u* Ai • u 


Ian 


Using Proposition 4.12, we have 


u* = 0 in fl, 


since u* is a solution to (4.20). 

Actually, we can see that there is no need for an extraction, since 0 is the only 
accumulation point for with respect to the weak topology. If we consider a 
subsequence it is still bounded in and therefore, using the same 

argument as above, zero is an accumulation point of this subsequence. For the 
strong convergence, we use (|4.25) to get 


[ [ vC/(^)-VU*+ [ 

Jn Jn Jan 


Ai • u. 


(4.26) 


Since 0, the right-hand side of (4.26) goes to zero when rj goes to zero. 

Hence, 

11 ^ 1 ( 0 )—>0 as 77 ^ 0 . 

□ 


Now, we take to be the solution of (4.23) and define the approximated 
resistivity (inverse of the conductivity) by 


(Jti 


Since 


Proposition 


1 

0-* 


+ Ai 


|VU* + Ai 


(4.27) 


4.13 


shows that — is a good approximation for — in the L^-sense. 


(Jri 


O'* 
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Proposition 4.14. Let u* he the true conductivity and let an he defined hy (4-21} 
We have 


1 

1 

an 

< 7 * 


0 as r] ^ 0. 


L2(o) 


5 Numerical illustrations 

We set 0 = I(x,y) G ^ < l|. We take a conductivity a G as 


represented on Figure 5.1 The potential Ai is chosen as 


Ai(x) = 10-2(| + 1;-| + i) 
so that Bi is constant in space. 

5.1 Optimal control 


We use the algorithm presented in section 4.2.1 We set a step size equal to 8 • 10“’^ 
and (T(o) = 3 as an initial guess. After 50 iterations, we get the reconstruction 
shown in Figure 5.2 The general shape of the conductivity is recovered but the 


conductivity contrast is not recovered. Moreover, the convergence is quite slow. It 
is worth mentioning that using two nonparallel electric current densities does not 
improve significantly the quality of the reconstruction. 

5.2 Fixed-point method 


We use the algorithm described in section 4.2.2, but slightly modified. The operator 
Q defined by 

{a^V[a] + aAi) ■ 3 


G[cr] ■■= a- 


is replaced by 


G[a] := 


(VFH + Ai) • J 


IVFH+Ail" 

which is analytically the same but numerically is more stable. Since the term 
|VF[cj] + Aip can be small, we smooth out the reconstructed conductivity a^n) 
at each step by convolving it with a Gaussian kernel. This makes the algorithm less 


unstable. The result after 9 iterations is shown in Figure 5.3 The convergence is 
faster than the gradient descent, but the algorithm still fails at recovering the exact 
values of the true conductivity. 
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Figure 5.1: Conductivity to be reconstructed. 
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Figure 5.2: Conductivity reconstructed by the optimal control method. 
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Figure 5.3: Conductivity reconstructed by the fixed point method. 
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5.3 Orthogonal field method 


We set r/ = 5 • 10 ^ and perform the computation described in section 4.2.3 The 


result we get is shown in Figure [5^ It is a scaled version of the true conductivity cr*, 
which means that the contrast is recovered. So assuming we know the conductivity 
in a small region of (or near the boundary dfl) we can re-scale the result, as 


shown in Figure 5.5 When tj goes to zero, the solution of (4.23) converges to the 


true potential 14 up to a scaling factor which goes to inhnity. When tj is large, the 
scaling factor goes to one but the solution becomes a ’’smoothed out” version 
of 14- This method allows an accurate reconstruction of the conductivity by solving 
only one partial differential equation. It covers the contrast accurately, provided we 
have a little bit of a prior information on ct*. 

Finally, we study the numerical stability with respect to measurement noise of 
the orthogonal field method. We compute the relative error defined by 


e := 


\^ri 


averaged over 150 different realizations of measurement noise on J. The results are 


shown in Figure 5.6 We show the results of a reconstruction with noise level of 2% 
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Figure 5.5: Conductivity recovered by the orthogonal field method after scaling. 
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Figure 5.6: Relative error with respect to measurement noise. 


(resp. 10%) in Figure [5?7| (resp. FigureClearly, the orthogonal method is quite 
robust with respect to measurement noise. 

6 Concluding remarks 

In this paper we have presented a new mathematical and numerical framework for 
conductivity imaging using magnetoacoustic tomography with magnetic induction. 
We developed three different algorithms for conductivity imaging from boundary 
measurements of the Lorentz force induced tissue vibration. We proved convergence 
and stability properties of the three algorithms and compared their performance. 
The orthogonal field method performs much better than the optimization scheme 
and the fixed-point method in terms of both computational time and accuracy. In¬ 
deed, it is robust with respect to measurement noise. In a forthcoming work, we 
intend to generalize our approach for imaging anisotropic conductivities by magne¬ 
toacoustic tomography with magnetic induction. 
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Figure 5.7: Reconstruction with the orthogonal field method with measurement 
noise level of 2%. 
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Figure 5.8: Reconstruction with the orthogonal field method with measurement 
noise level of 10%. 
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